Data are available in two formats:

1. Fishing effort by flag state and gear type at 100th degree resolution (i.e. too large to upload)
2. Fishing effort by MMSI at 10th degree resolution (i.e. too large to upload)

Each version of the fishing effort dataset is accompanied by a table of vessel information (e.g. gear type, flag state, dimensions, etc. listed in the metadata folder).

Note: When attempting to download multiple files, you will receive an email containing download links after the data has been packaged. For large requests, it may take up to 24 hours to receive this email.

Versions

The fishing effort data is updated periodically to include more recent data and as improvements are made to our models and input data. Data from different major/minor versions should not be combined, as each version uses a different set of inputs.

Current version

2.0: Fishing effort data for 2012-2020 using the latest AIS algorithms, neural network models, and vessel registry database.

Load Packages

We will use several R packages to work with the data. The tidyverse package contains several helpful packages for data manipulation (dplyr, tidyr, readr) and plotting (ggplot2), lubridate() which is great for working with dates. The tidymodels() package is key to this analysis for setting up all the models for machine learning.

Some of the Models used include :

PCA Lasso and Elastic-Net Regularized Generalized Linear Decision Tree Random Forest Boosted Tree KNN Neural Network

Shapefiles

Shapefiles of land and Exclusive Economic Zone (EEZ) boundaries are useful for providing context to maps of fishing efforts (too large to upload). For this analysis, we’ll use a land shapefile from the maps R package and an EEZ shapefile from MarineRegions.org. We’ll also load/convert both of these shapefiles to sf objects. The sf R package is a great package for working with vector data stored in dataframes, providing numerous functions for performing spatial operations.

World polygons from the maps package

# load("world.Rda")

Load EEZ polygons

# load("eezs.Rda")

Principal components analysis (PCA)

Principal components analysis is an ordination method allowing us to glean as much about our multivariate data as possible in a simplified number of dimensions. Here, we’ll use our dataset of fishing vessels to explore variable relationships and clustering by geartype in a PCA biplot. For this example with PCA, we will only use numeric predictors (i.e.mmsi,vessel_class_inferred_score,length_m_inferred,engine_power_kw)*.

We use the R built-in function prcomp() to calculate the principal components of the dataset. It’s important to specify scale = TRUE so that each of the variables in the dataset are scaled to have a mean of 0 and a standard deviation of 1 before calculating the principal components.

Also note that eigenvectors in R point in the negative direction by default, so we’ll multiply by -1 to reverse the signs.

# install.packages("imputeTS")

library(imputeTS) # for na.replace
#calculate principal components
vessels_pca <- vessels_sub %>%
  select_if(is.numeric) %>% 
  imputeTS::na.replace(0) %>% 
  scale() %>% 
  prcomp()



#reverse the signs
vessels_pca$rotation <- -1*vessels_pca$rotation

#display principal components
vessels_pca$rotation

#reverse the signs of the scores
vessels_pca$x <- -1*vessels_pca$x

#display the first six scores
# head(vessels_pca$x)

library(pca3d) # 2d & 3d pca plots
library(viridis) #color palette
library(grDevices)
# Variance explained by each PC
screeplot(vessels_pca, type = "lines")

screeplot(vessels_pca, type = "barplot")

pca2d(
   vessels_pca,
  components = 1:2,
  col = NULL,
  title = NULL,
  new = FALSE,
  axes.color = "black",
  bg = "white",
  radius = 1,
  shape = 18,
  palette = palette(hcl.colors(3, "viridis")) )

pca_3d <- pca3d(
  vessels_pca,
  components = 1:3)

## [1] 0.7222141 0.5503004 0.3397774
## Creating new device

## a bit more fancy:
  ## black background, white axes,
  ## centroids
pca2_3d <- pca3d(
  vessels_pca,
    fancy= TRUE, bg= "black", 
    axes.color= "white", new= TRUE) 

## [1] 0.7222141 0.5503004 0.3397774
## Creating new device

The PCA analysis indicates that PC1 & PC2 represent a good amount of the data compared to the rest.

Fig3. 3D image of the firts 3 PCA’s

Fig4. 3D image of the firts 3 PCA’s

Data splitting

set.seed(3435)

# Split the dataset into train and test

vessels_split <- initial_split(vessels_sub, prop = 0.8, strata = vessel_class_gfw)
vessels_train <- training(vessels_split)
vessels_test <- testing(vessels_split)
# dim(vessels_train)

In order to have a good training set, we set up a prop equal to 80% and our training set has 29 variables with 7998 obs.

Training set Visualisation

1. Corrplot

vessels_train %>% 
  select_if(is.numeric ) %>% 
  dplyr::select(- vessel_class_inferred_score) %>% 
  dplyr::select(- starts_with("fishing")) %>% 
  cor(use = "complete.obs") %>% 
  corrplot(type = "lower", diag = FALSE)

From the corrplot, we can see that tonnage_gt_ registry and engine_power_kw_gfw are positively correlated but length_mregitery is negatively correlated with engine power.

2. Barplot

vessels_train %>% 
  ggplot(aes(x = vessel_class_gfw)) +
  geom_bar(aes(fill=vessel_class_gfw))

From the barplot of the chosen predictors : “set_longlines”, “set_gillnets”, “tuna_purse_seines”, “fishing”,
“other_purse_seines”,” drifting_longlines”. We notice that fishng class has more observations than the others.

Cross-Validation folds

For the cross-validation, I’m using v = 3 folds and stratifying the outcome variable for geartype vessel_class_gfw.

set.seed(3435)
vessels_folds <- vfold_cv(data = vessels_train, v = 3, strata = vessel_class_gfw)

Using a subset of the data and given that the training set only has 29 vars and many observations. we chose to use lower v-folds to avoid having an imbalance in a given class where an outcome variable may either be present or not.

Elastic Tuning Net

We’ll be fitting and tuning an elastic net, tuning penalty and mixture using multinom_reg with the glmnet enginefor classification. We set up a model and workflow, create a regular grid for penalty and mixture with 10 levels each; mixture should range from 0 to 1. The penalty range used is from -5 to 5 (it’s log-scaled).

# Elastic Net Model
elastic_net_spec <- multinom_reg(penalty = tune(), 
                                 mixture = tune()) %>% 
  set_mode("classification") %>% 
  set_engine("glmnet")

en_workflow <- workflow() %>% 
  add_recipe(vessels_recipe) %>% 
  add_model(elastic_net_spec)

en_grid <- grid_regular(penalty(range = c(-5, 5)), 
                        mixture(range = c(0, 1)), levels = 10)

Using library(tictoc) to print the system running time for my model

library(tictoc)
tic()
tune_res_1 <- tune_grid(
  en_workflow,
  resamples = vessels_folds, 
  grid = en_grid
)

# save(tune_res_1, file ="tune_res.Rda")
# load("tune_res.Rda")
ggplot2::autoplot(tune_res_1)

toc()

## 60.752 sec elapsed

We’re fitting a grid of 100 combinations of values for the hyperparameters penalty and mixture, or 100 models, and we’re fitting each of these once per fold, for a total of 200 models.

tune_res_1 <- tune_grid(
  en_workflow,
  resamples = vessels_folds, 
  grid = en_grid
)

autoplot(tune_res_1)

Smaller values of penalty, or “amount of regularization,” and smaller values of mixture, or “proportion of lasso penalty,” tend to result in higher ROC-AUC and accuracy values.

Model Prediction, fitting & Visualization

best_model <- select_best(tune_res_1, metric = "roc_auc")

en_final <- finalize_workflow(en_workflow, best_model)

en_final_fit <- fit(en_final, data = vessels_train)

predicted_data <- augment(en_final_fit, new_data = vessels_test) %>% 
  dplyr::select(vessel_class_gfw, starts_with(".pred")) 


predicted_data %>% 
  conf_mat(truth =vessel_class_gfw, estimate = .pred_class) %>%
  autoplot(type = "heatmap")

In the heat map, the correct classifications are always on the diagonal and misclassifications on the off-diagonals. The model correctly identified fishing geartype, for example. From the heat map, we can see that the model doesn’t seem to be able to distinguish other_purse_seines well at all; it misclassifies other geartypes as fishing frequently.

predicted_data %>% roc_auc(vessel_class_gfw, .pred_fishing:.pred_tuna_purse_seines)

## # A tibble: 1 × 3
##   .metric .estimator .estimate
##   <chr>   <chr>          <dbl>
## 1 roc_auc hand_till      0.757

roc_auc in this model is around 72%.

Decision Trees

One of our favorite machine learning algorithms is decision trees. A single decision tree generally doesn’t have great “out of the box” model performance, and even with considerable model tuning, they are unlikely to perform as well as other approaches. Decision trees are non-parametric, meaning they do not make any assumptions about the underlying data-generating process.

At their core, decision trees work by splitting the features into a series of yes/no decisions. These splits divide the feature space into a series of non-overlapping regions, where the cases are similar in each region. To understand how a given prediction is made, one simply “follows” the splits of the tree (a branch) to the terminal node (the final prediction). This splitting continues until a specified criterion is met.

Decision tree hyperparamters

Decision trees have three hyperparameters as shown below. These are standard hyperparameters and are implemented in {rpart}, the engine used to fit decision tree models in the previous section. Alternative implementations may have slightly different hyperparameters (see the documentation for parsnip::decision_tree() details on other engines).

First,I set up a decision tree model and workflow. Tune the cost_complexity hyperparameter using 10 levels, range = c(-3, -1) and specify that the metric we want to optimize is roc_auc. Lastly, print an autoplot() of the results.

tree_spec <- decision_tree() %>%
  set_engine("rpart")

class_tree_spec <- tree_spec %>%
  set_mode("classification")



class_tree_wf <- workflow() %>%
  add_recipe(vessels_recipe) %>%  
  add_model(class_tree_spec %>% set_args(cost_complexity = tune())) 

param_grid_tr <- grid_regular(cost_complexity(range = c(-3, -1)), levels = 10)



tune_res_tr <- tune_grid(
  class_tree_wf, 
  resamples = vessels_folds, 
  grid = param_grid_tr, 
  metrics = metric_set(roc_auc)
)


autoplot(tune_res_tr)

A single decision tree performs better with larger complexity penalty( peak at mid level complexity)

tree_metrics <- collect_metrics(tune_res_tr) %>% 
arrange(desc((mean)))

best_complexity <- select_best(tune_res_tr )

Decision Tree Visualization

Using rpart.plot, I fit and visualize my best-performing pruned decision tree with the training set.

class_tree_final <- finalize_workflow(class_tree_wf, best_complexity)

class_tree_final_fit <- fit(class_tree_final, data = vessels_train)

class_tree_final_fit %>%
  extract_fit_engine() %>%
  rpart.plot()

Random forest

In this section, I set up a random forest model with rand_forest(), a workflow using the ranger engine and set importance = "impurity". Then, once again I will tune mtry, trees, and min_n.

Random forest tree hyperparameters

I start by creating a regular grid with 8 levels each. For each hyperparameter (i.e. this model has 3 tuning parameters I specify the following ranges :

mtry: # Randomly Selected Predictors reduce variance because mtry allows to choose vars randomly. Depends on the number of columns. mtry = P/3 or sqrt P Note that mtry should not be smaller than 1 or larger than 8.

trees: # Trees (type: integer, default: 500L)

min_n: (stopping criterion) defines the minimal node size (node 1 and has 10 obs). min_n depends on the mode. For regression, a value of 5 is the default. For classification, a value of 10 is used.

class_rf_model = rand_forest(mtry = tune(), trees = tune(), min_n = tune()) %>% 
  set_engine("ranger", importance = "impurity") %>% 
  set_mode("classification") 


class_rf_wflw = workflow() %>% 
  add_recipe(vessels_recipe) %>% 
  add_model(class_rf_model) 


library(tictoc)
tic()

rf_grid <- grid_regular(
  mtry(range = c(1, 8)), #set mtry to 1-8
  trees(range = c(100,200)), #set trees 100 to 200
  min_n(range = c(2, 10)), # set min_n between 2 and 20
  levels = 8) #what's a good level?

rf_tune <- tune_grid( 
  class_rf_wflw,
  resamples = vessels_folds, 
  grid = rf_grid, 
  metrics = metric_set(roc_auc) 
  )
toc()

## 4748.907 sec elapsed

autoplot(rf_tune)

saveRDS(rf_tune,file = "rf_tune.RDS")
rf_tune_results <- readRDS("rf_tune.RDS")

rf_tree_metrics <- collect_metrics(rf_tune) %>% 
arrange(desc((mean)))
rf_tree_metrics

## # A tibble: 512 × 9
##     mtry trees min_n .metric .estimator  mean     n std_err .config             
##    <int> <int> <int> <chr>   <chr>      <dbl> <int>   <dbl> <chr>               
##  1     6   200    10 roc_auc hand_till  0.827     3 0.00502 Preprocessor1_Model…
##  2     6   185    10 roc_auc hand_till  0.827     3 0.00385 Preprocessor1_Model…
##  3     7   200    10 roc_auc hand_till  0.826     3 0.00458 Preprocessor1_Model…
##  4     8   185    10 roc_auc hand_till  0.826     3 0.00358 Preprocessor1_Model…
##  5     5   185    10 roc_auc hand_till  0.826     3 0.00395 Preprocessor1_Model…
##  6     6   171    10 roc_auc hand_till  0.826     3 0.00472 Preprocessor1_Model…
##  7     7   171    10 roc_auc hand_till  0.826     3 0.00467 Preprocessor1_Model…
##  8     7   185     8 roc_auc hand_till  0.826     3 0.00454 Preprocessor1_Model…
##  9     6   200     7 roc_auc hand_till  0.826     3 0.00446 Preprocessor1_Model…
## 10     5   200    10 roc_auc hand_till  0.826     3 0.00466 Preprocessor1_Model…
## # … with 502 more rows

best_complexity_rf <- select_best(rf_tune)

Random forest Predictions

Using vip(), with your best-performing random forest model fit on the training set.

Which variables were most useful? Which were least useful? Are these results what you expected, or not?

From the above graph, special attack is the most important variable. Attack, speed,defense, hp,and special are equally as important.

rf_final <- finalize_workflow(class_rf_wflw, select_best(rf_tune)) 

rf_fit <- fit(rf_final, vessels_train)  

 vip_rf <- rf_fit %>% 
  extract_fit_engine() %>%
  vip()
 
 vip_rf

Boosted tree

Here, I set up a boosted tree model and workflow using the xgboost engine. Then, I use tune() to set a boosted model trees. The grid used below contains 10 levels; where trees range from 10 to 2000. Lastly, I specify roc_auc in my model metrics and print an autoplot() to view the results.

boost_model = boost_tree(trees = tune()) %>% #set up boosted tree model
  set_engine("xgboost") %>% # xgboost engine calculates lambda #high lambda is good 
  set_mode("classification") 


boost_wflw = workflow() %>%
  add_recipe(vessels_recipe) %>% 
  add_model(boost_model) 

boost_grid <- grid_regular(
  trees(range = c(10,2000)), 
  levels = 10) 

boost_tune_res <- tune_grid( 
  boost_wflw, 
  resamples = vessels_folds, 
  grid = boost_grid, 
  metrics = metric_set(roc_auc) 
  )

# saveRDS(boost_tune_res,file = "boost_tune.RDS")
# boost_results <- readRDS("boost_tune.RDS")
# autoplot(boost_results)

autoplot(boost_tune_res) 

ROC AUC peaks 231 then slowly declines again.

boost_metrics <- collect_metrics(boost_tune_res) %>% 
arrange(desc((mean)))

best_complexity_boost <- select_best(boost_tune_res)

The roc_auc of the best-performing boosted tree model on the folds

Decion Tree, Random Forest, & Boosting model predictions

models = full_join(tree_metrics, rf_tree_metrics)
joined_models = full_join(models, boost_metrics)
joined_models 

## # A tibble: 532 × 10
##    cost_complexity .metric .estimator  mean     n std_err .config     mtry trees
##              <dbl> <chr>   <chr>      <dbl> <int>   <dbl> <chr>      <int> <int>
##  1         0.001   roc_auc hand_till  0.776     3 0.00610 Preproces…    NA    NA
##  2         0.00167 roc_auc hand_till  0.771     3 0.00679 Preproces…    NA    NA
##  3         0.00278 roc_auc hand_till  0.753     3 0.0130  Preproces…    NA    NA
##  4         0.00464 roc_auc hand_till  0.737     3 0.00481 Preproces…    NA    NA
##  5         0.00774 roc_auc hand_till  0.728     3 0.00110 Preproces…    NA    NA
##  6         0.0129  roc_auc hand_till  0.725     3 0.00319 Preproces…    NA    NA
##  7         0.0215  roc_auc hand_till  0.623     3 0.0225  Preproces…    NA    NA
##  8         0.0359  roc_auc hand_till  0.596     3 0.00437 Preproces…    NA    NA
##  9         0.0599  roc_auc hand_till  0.596     3 0.00437 Preproces…    NA    NA
## 10         0.1     roc_auc hand_till  0.596     3 0.00437 Preproces…    NA    NA
## # … with 522 more rows, and 1 more variable: min_n <int>

The best model is the random forest model, with an AUC ROC of 0.628.

Best-performing model on the testing set

rf_test_fit = augment(rf_fit, new_data = vessels_test) 
roc_auc(rf_test_fit, truth = vessel_class_gfw, .pred_fishing:.pred_tuna_purse_seines) 

## # A tibble: 1 × 3
##   .metric .estimator .estimate
##   <chr>   <chr>          <dbl>
## 1 roc_auc hand_till      0.818

AUC value of the best-performing model on the testing set.

autoplot(roc_curve(rf_test_fit, truth = vessel_class_gfw, .pred_fishing:.pred_tuna_purse_seines)) 

conf_mat(rf_test_fit, truth = vessel_class_gfw, estimate = .pred_class) %>% 
  autoplot(type = "heatmap") + 
  theme(axis.text.x = element_text(angle = 90, hjust=1)) 

K-nearest neighbor

library(kknn)
knn_spec <- 
  nearest_neighbor(neighbors = 4) %>% # we can adjust the number of neighbors 
  set_engine("kknn") %>% 
  set_mode("classification") 

knn_wflow <-
 workflow() %>%
 add_recipe(vessels_recipe) %>% 
 add_model(knn_spec)

knn_res <- 
  knn_wflow %>% 
  fit_resamples(
    resamples = vessels_folds, 
    metrics = metric_set(
      recall, precision, f_meas, 
      accuracy, kap,
      roc_auc, sens, spec),
    control = control_resamples(save_pred = TRUE)
    ) 

knn_res %>% collect_metrics(summarize = TRUE)

## # A tibble: 8 × 6
##   .metric   .estimator  mean     n std_err .config             
##   <chr>     <chr>      <dbl> <int>   <dbl> <chr>               
## 1 accuracy  multiclass 0.423     3 0.00554 Preprocessor1_Model1
## 2 f_meas    macro      0.423     3 0.00367 Preprocessor1_Model1
## 3 kap       multiclass 0.196     3 0.00814 Preprocessor1_Model1
## 4 precision macro      0.422     3 0.00863 Preprocessor1_Model1
## 5 recall    macro      0.426     3 0.00931 Preprocessor1_Model1
## 6 roc_auc   hand_till  0.729     3 0.00445 Preprocessor1_Model1
## 7 sens      macro      0.426     3 0.00931 Preprocessor1_Model1
## 8 spec      macro      0.838     3 0.00188 Preprocessor1_Model1

# knn_results <- saveRDS(boost_tune_res,file = "boost_tune.rds")
# knn_results <- readRDS("boost_results.rds")

The roc_auc of the best-performing knn here is 0.40

Neural network

#  # ?keras::install_keras()
# # install_keras(Tensorflow = "1.13.1",
# # restart_session = FALSE)
# library(keras)
# library(tensorflow)
# 
# nnet_spec <-
#   mlp() %>%
#   set_mode("classification") %>% 
#   set_engine("keras", verbose = 0) 
# 
# nnet_wflow <-
#  workflow() %>%
#  add_recipe(vessels_recipe) %>% 
#  add_model(nnet_spec)
# 
# nnet_res <- 
#   nnet_wflow %>% 
#   fit_resamples(
#     resamples = vessels_folds, 
#     metrics = metric_set(
#       recall, precision, f_meas, 
#       accuracy, kap,
#       roc_auc, sens, spec),
#     control = control_resamples(save_pred = TRUE)
#     )